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ABSTRACT 

A Monte Carlo code to solve the transfer of Lya photons is developed, which can 
predict the Lya image and two-dimensional Lya spectra of a hydrogen cloud with 
any given geometry, Lya emissivity, neutral hydrogen density distribution, and bulk 
velocity field. We apply the code to several simple cases of a uniform cloud to show 
how the Lya image and emitted line spectrum are affected by the column density, 
internal velocity gradients, and emissivity distribution. We then apply the code to two 
models for damped Lya absorption systems: a spherical, static, isothermal cloud, and a 
flattened, axially symmetric, rotating cloud. If the emission is due to fluorescence of the 
external background radiation, the Lya image should have a core corresponding to the 
region where hydrogen is self-shielded. The emission line profile has the characteristic 
double peak with a deep central trough. We show how rotation of the cloud causes 
the two peaks to shift in wavelength as the slit is perpendicular to the rotation axis, 
and how the relative amplitude of the two peaks is changed. In reality, damped Lya 
systems are likely to have a clumpy gas distribution with turbulent velocity fields, which 
should smooth the line emission profile, but should still leave the rotation signature of 
the wavelength shift across the system. 

Subject headings: line:formation - radiative transfer - scattering - quasars: absorption 
lines 



1. Introduction 

The high column density absorption systems (damped Lya and Lyman limit systems) are a key 
part of understanding galaxy formation. As a galaxy collapses from the highly ionized intergalactic 
medium, the gas will inevitably go through a phase of a self-shielded cloud of atomic hydrogen before 
it can cool and collapse further into molecular clouds, form stars, and increase its metallicity. Many 
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of the damped Lya systems (DLAs) at high redshift may be gaseous halos in the process of forming 
galaxies rather than fully formed spiral disks, especially in view of the very low metallicities at high 
redshifts (Lu et al. 1996; Prochaska, Gawiser, k Wolfe 2001). If so, the internal structure of these 
objects should reveal to us the detailed processes by which galaxies form. Several basic questions 
arise in relation to this: how big are the absorption systems? Is the gas smoothly distributed, or is 
it in the form of clumps and an interclump medium? If so, how big are the clumps? What is the 
dynamical state of the gas? How does its mean rotation velocity compare to its velocity dispersion? 

Some of these questions may be addressed by studying the associated metal absorption lines 
(e.g., Prochaska k Wolfe 1997, 1998) or by radio and optical observations of galaxies found to be 
associated with the absorbers (see Briggs et al. 1989; Djorgovski et al. 1996), although the majority 
of DLAs at high redshifts might not be associated with luminous galaxies. 

An alternative observational probe of the structure of DLAs and Lyman limit systems (LLSs) 
may be found in their Lya emission produced in hydrogen recombinations and the subsequent 
scattering of Lya photons through the gas (Hogan k Weymann 1987; Gould k Weinberg 1996). 
The source of the ionization may be the external cosmic ionizing radiation, shock-heating due to 
the gravitational collapse of the gas in a dark matter halo, or star formation within the system. In 
the case of fluorescence of external radiation, the maximum surface brightness is achieved in any 
system with AT HI > 10 18 cm -2 , where all the incident external photons are absorbed. 

Efforts have been made to image known DLAs and LSSs towards high-z quasars. DLAs 
towards PKS0528-250 (Warren k M0ller 1996) and Q0151+048A (Fynbo, M0ller, k Warren 1999) 
and one LLS towards Q1205-30 (Fynbo, Thomsen, k Warren 2000) have been successfully detected 
using tuned narrow band filters. There are also some spectroscopically confirmed detections of Lya 
emission from DLAs (e.g. Hunstead, Pettini, k Fletcher 1990; Djorgovski 1998; Pettini et al. 1995; 
Leibundgut k Robertson 1998). Most of these absorber systems with Lya emission detected are 
near the redshift of the background quasar, which suggests a Lya source due to photoionization by 
the quasar (Fynbo, M0ller, k Warren 1999; Fynbo, Thomsen, k Warren 2000). 

Hydrogen Lya line is a resonance line. The problem of radiative transfer of resonance line 
radiation can be approximately solved analytically under certain conditions (e.g., Harrington 1973, 
1974; Neufeld 1990). These analytic solutions can be found only for a limited number of cases such 
as a static, extremely opaque and plane-parallel medium. In some cases, numerical methods are 
used to solve the transfer equation (e.g., Adams 1972; Hummer k Kunasz 1980; Urbaniak k Wolfe 
1981). For a more general geometry, density distribution, and kinematics, the Monte Carlo method 
turns out to be very useful and thus it has been applied to many problems (e.g., Auer 1968; Avery 
k House 1968; Caroff, Noerdlinger, k Scargle 1972; Panagia k Ranieri 1973; Bonilha et al. 1979; 
Natta k Beckwith 1986; Ahn, Lee, k Lee 2000, 2001, 2002). 

This paper presents the results of a numerical Monte Carlo code we have developed to compute 
the transfer of Lya photons through an arbitrary distribution of hydrogen. The code predicts the 
two-dimensional (2-D) image and line-spectrum along any given observed direction, for a general 
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three-dimensional distribution of gas of given geometry, Lya emissivity, neutral hydrogen density, 
and bulk velocity field. In this paper, it will be applied only to spherical and axially symmetric 
clouds, although in the future it can be applied to results of numerical simulations of galaxy 
formation. 

We describe the Monte Carlo code in §2. In §3, we apply the code to several simple cases of 
hydrogen clouds to demonstrate the effects of the neutral hydrogen column density, Lya emissivity 
and bulk velocity field on the spectra of escaped Lya photons. In §4, we model DLAs as static and 
rotating clouds and investigate their 2-D Lya emission spectra. Finally, we have a brief summary 
and discussion in §5. 

2. Monte Carlo Simulation of Lya Scattering 

2.1. Description of the Method 

The method used to numerically compute the images and spectra of Lya emission consists 
of generating random realizations of the trajectory of a large number of Lya photons as they 
are scattered within the specified gas distribution. In a nutshell, the method proceeds through 
the following steps: First, we generate the point where the photon is created, with the emissivity 
distribution proportional to the recombination rate at every point. Second, the optical depth r that 
the photon will travel through before it is scattered is randomly generated with a probability density 
e~ r . The spatial location of the scattering at optical depth r from the point of emission is then 
determined along a randomly chosen direction, with the knowledge of the neutral hydrogen density 
distribution and the scattering cross section (see eq.[2] below). We then generate the velocity of 
the atom that scatters the photon, as described below, compute the new frequency of the photon 
after scattering, and generate the new direction of the photon. The process of propagation and 
scattering is repeated until the photon escapes the modeled system. 

As these random realizations of photon trajectories are being carried out, a Lya image of the 
system along a specified direction can be created. The image consists of a three-dimensional array 
of the two projected coordinates perpendicular to the direction of the image and the frequency of 
the escaped photons. The array contains the mean number of photons emitted at every projected 
position and frequency. At every scattering of a photon, the probability that the photon is re- 
emitted in the direction of the image and escapes is separately calculated. This probability is 
added on the element of the image array corresponding to the projected position and frequency 
of the photon. In practice, for the majority of the scatterings the optical depth is large and the 
probability of escape is negligible, so one can avoid the computation of the optical depth in the 
direction of the image on most scatterings. 

We now describe how the optical depth of the Lya photons is computed and the velocity of 
the atom at each scattering is generated. The scattering cross section of Lya photons as a function 
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of the frequency in the rest frame of the hydrogen atom is 

7re 2 Av L /2n 

where /12 = 0.4162 is the Lya oscillator strength, uq = 2.466 x 10 15 Hz is the line center frequency, 
Aul = 4.03 x 10~ 8 z;o = 9.936 x 10 7 Hz is the natural line width, and the other symbols have their 
usual meaning. For a hydrogen atom with a velocity component v z along the photon's direction in a 
fixed "laboratory" frame, the term v — uq in the above equation becomes (i/j — vq) — {v z /c)vq, where 
Vi is the frequency of the incident photon in the "laboratory" frame. For a Maxwellian distribution 
of atom velocities, the resulting average cross section is 

<r(si) = /i2 H(a, Xi ) , (2) 

m e cAuo 



where 



a [ + °° e~ y2 



is the Voigt function, Avp = {v p /c)vq is the Doppler frequency width, v p = (2/cT/m//) 1 / 2 is the 
atom velocity dispersion times \/2, T is the gas temperature, inn is the mass of the hydrogen 
atom, Xi = (i/j — vq)/ Avd is the relative frequency of the incident photon in the laboratory frame, 
and a = Avl/ \2Aud) is the relative line width. The optical depth along the photon trajectory is 
computed by integrating the cross section in equation (2) times the atomic hydrogen density. 

Once a spatial location where the photon is scattered has been determined, the velocity v z 
along the direction of the incident photon of the atom responsible for the scattering obeys the 
following distribution: 

f( u z) = -- ( n' 2 ■ 2 H-\a,Xi) , (4) 

7T (Xi - u z ) 1 + a 1 

where u z = v z /v p . We describe how we generate random numbers with this distribution in the 
Appendix. The two velocity components perpendicular to the direction of the photon simply follow 
a Gaussian distribution parameterized by the local temperature. 

After the total velocity of the atom is chosen according to the above distribution, we first 
perform a Lorentz transform of the direction and frequency of the photon to the rest frame of the 
atom. The direction of the scattered photon is then generated according to a dipole distribution 
in this frame. Its frequency differs from the incident one by the recoil effect which is taken into 
account in the code (it is negligibly small in the applications we will present in this paper). The 
direction and frequency of the scattered photon is then transformed back to the laboratory frame, 
a new optical depth is chosen, and the entire process is repeated until the photon escapes from 
the cloud. In the presence of a fluid velocity, the procedure changes only by replacing the incident 
frequency Xi in equations (2) and (4) by its value in the fluid frame, xa = x% — [vi z / o){vq / Au£>) , 
where Vf z is the fluid velocity parallel to the photon's direction (valid in the non-relativistic regime). 
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As mentioned previously, the Lya image and spectrum of the system are calculated as the 
photon trajectories are randomly generated. At every photon scattering, we compute first the 
optical depth for the photon to escape the cloud along the direction of the image, and the photon 
is added to the corresponding element of the image array in projected position and frequency, with 
the weight e _T (l + fi 2 )dQ, where r is the optical depth for escaping, \i is the cosine of the angle 
between the incident photon and the image direction, and dO is the solid angle subtended by the 
pixel (both [i and dSl are measured in the rest frame of the atom). The factor 1 + fi 2 accounts 
for the dipole probability distribution of the direction of the photon after scattering. A photon 
also has a probability to directly escape the system along the direction of the image when it is 
created which is added to the corresponding element of the image array. After the contribution to 
the image has been added, a random direction of the new scattered photon is generated according 
to the same dipole distribution, and the realization of the photon trajectory is continued until the 
photon finally escapes from the cloud. The total number of photons in the simulation is chosen to 
give an acceptable signal-to-noise ratio (the typical number of photons used for simulations in this 
paper is on the order of ~ 10 3 — 10 ). 



2.2. Tests of the Code 




Fig. 1. — Comparison between results from our Monte Carlo simulations and the analytic solutions. 
Dotted lines are the analytic solutions (Neufeld 1990) of the Lya spectra at one boundary of a 
slab with a midplane source and different scattering optical depths. Solid lines are those from 
simulations. In the left panel, the recoil effect is neglected in the simulations; in the right panel, 
this effect is taken into account in the simulations (see the text). 
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We test our numerical code for the case of a static, plane-parallel slab, for which Neufeld (1990) 
derived an analytic solution of the mean intensity in the limit of a large optical depth. We simulate 
the case of a midplane source radiating line-center Lya photons in a plane-parallel slab of uniform 
density for three optical depths: tq = 10 4 , 10 5 , and 10 6 , where tq = cr(xi = 0)Ahi is the line-center 
optical depth from the midplane to the boundary of the slab. We assume a temperature T = 10K 
(a « 1.49 x 1CT 2 ). As we have mentioned, the recoil effect is taken into account in the code. For 
the purpose of comparison, we first turn off this effect and run the simulation. We compare the 
spectra of the escaped Lya photons with the analytic results (eq.[2.24] in Neufeld 1990; note that 
his definition of r differs from ours by a factor of y/ir as a <C 1. See also Ahn, Lee, & Lee 2001). 
The left panel of Figure 1 shows excellent agreement, becoming better as the optical depth increases 
as expected since the analytic solution applies in the limit of a large optical depth (^/tttq > 10 3 /a, 
Neufeld 1990). We then turn on the recoil effect and rerun the simulation. The results are shown 
in the right panel of Figure 1. The double-peaked frequency distribution of the escaped photons 
becomes asymmetric - more photons appear at the red part. The profile from analytic solutions 
can be regarded as an average of the red part and the blue part. The magnitude of the recoil 
effect is easily understood as reflecting the thermalization of photons around frequency v inside 
the cloud (Wouthuysen 1952; Field 1959), which modifies their abundance by a factor of e~ x / XT , 
where x = Au/Au D and xt = kT/(hAu D ) = (m H c 2 kT '/2) 1 / 2 / '(hu ). Obviously, the recoil effect 
becomes more important for low temperature and large optical depth. 

We also test the code for a spherically-symmetric case with a particular bulk velocity field. 
Loeb & Rybicki (1999) study and simulate the scattered Lya radiation around a point source before 
cosmological reionization. The Lya photons are scattered by the intergalactic medium (IGM) which 
is undergoing Hubble expansion. According to their approximation, the gas temperature is regarded 
as zero. We run a corresponding simulation with our code and find that the frequency distribution 
and the surface brightness profile of the escaped Lya photons shown in Loeb & Rybicki (1999) are 
well reproduced. 

3. Spherical Clouds of Uniform Density 

We start applying our code to the most simple case of spherical hydrogen clouds with uniform 
density. We consider the following cases: 

Case 1: uniform emissivity, static cloud 

Case 2: uniform emissivity, expanding or contracting cloud 

Case 3: central point source, static cloud 

Case 4: central point source, expanding or contracting cloud 

For each case, we perform two runs with different column densities (2 x 10 18 cm~ 2 and 2 x 
10 20 cm~ 2 , corresponding to line-center optical depths tq = 8.3 x 10 4 and 8.3 x 10 6 , respectively). 
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Fig. 2. — Frequency distribution of Lya photons escaped from a uniform hydrogen cloud with 
neutral hydrogen column density 2 x 10 18 cm -2 for the cases of uniform emissivity and central 
point source. Different line types are for a static, expanding, or contracting cloud. 

The temperature is assumed to be 2 x 10 4 K everywhere. In the case of expansion, we set it to 
be Hubble-like (i.e., the velocity is proportional to the radius), with the velocity at the edge of the 
system fixed to be 200 km s -1 . The cases of contraction are done in the same way; their line spectra 
can generally be obtained by simply mirror-reflecting in frequency the photons escaped from the 
expanding cloud about the rest frame Lya frequency. 

Figures 2 and 3 show the distribution of escaped photons for all the cases. For static cases, as 
expected, the escaped photons have a double-peaked frequency distribution. Except for the very 
small effect of atomic recoil (see Field 1959), the two peaks are exactly symmetric; the differences 
in the figures are due to the simulation noise. 

The results obtained for the emergent spectrum can be understood by noting the typical tra- 
jectory in frequency and space followed by a photon (e.g., Osterbrock 1962; Adams 1972; Urbaniak 
& Wolfe 1981; Ahn, Lee, & Lee 2000). In the case of a Lyman limit system (like in Fig. 2), the pho- 
tons are always most likely to be scattered by atoms that have the right velocity along the photon 
direction so that, in their rest frame, the photon frequency is shifted very close to the line center. 
This implies there is practically no correlation between the frequency of the photon before and 
after scattering. At any random time, a photon is most likely to be found within a frequency range 
Az^d = {v p /c)vq of the line center, and will occasionally make larger excursions away from the line 
center when scattered by an atom having a large velocity perpendicular to the photon's direction. 
After any such excursion, the photon will likely return immediately to the line center with one or 
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Fig. 3. — Same as Fig. 2, but for neutral hydrogen column density 2 x 10 cm 



a few scatterings. Therefore, the photon escapes the cloud not by diffusing in frequency, but by 
making a single large jump when scattered by an atom in the high-velocity tail of the Maxwellian 
distribution, which reduces its optical depth for escaping to a value of order unity. For a Lyman 
limit system with Nm = 2 x 10 18 cm -2 , the line-center optical depth is tq = 8.3 x 10 4 , and the 
photons escape when their optical depth is r = T$e~ x2 ~ 1 (where x = Av/ Av^), which implies 
\x\ ~ v / hiTt) = 3.4. 

In the case of a DLA system, the scattering history differs from the Lyman limit case once 
the photon reaches a sufficiently large excursion in frequency to make it more likely that the next 
scattering is caused by a random atom far from the line center, rather than an atom moving with 
the right velocity to shift the photon frequency to the line center in the atom frame. Then, the 
evolution of the photon is described by diffusion in frequency. Since the photons are now reaching 
their escaping frequency not by a large jump but by a series of small steps, they undergo greater 
spatial diffusion than in Lyman limit systems. This, plus the fact that the optical depth has now 
a power-law instead of Gaussian dependence at large x, broadens the width of the two emission 
peaks. Of course, the emission peaks also move further from the line center because a smaller 
scattering cross section is required for the photons to escape. 

In Figure 4, we can see that the surface brightness produced by a central point source is more 
extended for a DLA system than for a Lyman limit system, owing to the greater spatial diffusion. 

For an expanding cloud, photons should escape on average with a redshift because they are 
doing work on the expansion of the cloud as they are scattered. For a Lyman limit system with 
a central source, a photon undergoing a large positive frequency jump will be moved back to the 
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Fig. 4. — Surface brightness profiles, in magnitudes per solid angle (with an arbitrary zero point). 
Curves are normalized in such a way that the total energy of the photons remains the same. In the 
right panel, the drop at r/R ~ 0.15 for the static case is due to simulation noise. 



line center relative to the fluid as it travels through the cloud, owing to the Hubble-like expansion 
of the cloud. On the other hand, a negative frequency jump will allow the photon to escape 
directly. In our case, the velocity at the cloud edge, 200 km s -1 , is much larger than the atomic 
velocity dispersion, so a photon needs to undergo many positive frequency jumps in order to diffuse 
spatially through the cloud and be able to escape on the blue side of the line, which explains why 
the blue peak is highly suppressed (see Fig. 3 of Ahn h Lee 1998 for a similar explanation). The 
situation is exactly reversed for a contracting cloud. The case of uniform emissivity broadens the 
line, essentially because of the different velocities of the emission sites of the photons; there is also 
the additional effect that photons emitted near the edge of the cloud are more likely to escape on 
the blue (red) peak for an expanding (contracting) cloud. In DLAs, the line is also broader for a 
point source compared to a Lyman limit system because of the power-law dependence of the cross 
section on frequency and the greater degree of spatial diffusion. 

To summarize, these simple models show how the frequency and spatial distributions of Lya 
photons escaped from a cloud are related to the Lya emissivity distribution, the bulk velocity field 
and the column density of the cloud. 
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4. Lya Emission from Model Damped Lya Systems 

We now apply our code to a gas cloud with an isothermal density profile as a model for DLAs. 
The nature of DLAs is still a subject of debate. They could be protogalaxies with a rotating disk 
component (e.g., Prochaska & Wolfe 1997, 1998) or spherically distributed clouds of gas moving 
randomly in halos (e.g., Mo 1994; Haehnelt, Steinmetz, & Rauch 1998; McDonald & Miralda- 
Escude 1999). The other mechanism that may give rise to DLAs is large-scale outflows due to 
galactic winds in dwarf galaxies (e.g., Nulsen, Barcons, & Fabian 1998; Tenorio-Tagle et al. 1999; 
Schaye 2001). In this paper, we focus on the first two pictures. We investigate the emergent spectra 
and spatial distribution of the Lya emission, to show what it can tell us about the internal structure 
of the system. 

4.1. Model Description 

We model the gas cloud producing a DLA assuming that it forms in a dark matter halo of 
mass lO n M0, virialized at redshift z = 3. In an Qm = 1 universe with Hq = 70kms~ 1 Mpc~ 1 , 
the corresponding virial radius is r v - a = 24kpc, and the virial velocity V w - U = 104 km s -1 (e.g., 
Padmanabhan 1993). We assume the gas density profile is a singular isothermal sphere with a 
cutoff at the virial radius and a fraction of the halo mass in gas of 0.05. 

We consider two different cases as an illustration to show how the rotation of a system could 
be probed by observations of Lya emission. The first case is a spherically symmetric cloud, and 
the second a rotating oblate ellipsoid with an axis ratio of 0.5 for the gas distribution. The 
rotating velocity is set to be V vot = ■\/2/3V v i I . In one case, we assume the gas to be static and 
at a temperature 2 x 10 4 K (a typical temperature of gas that cools after shock-heating but stays 
photoionized) . We also consider another case where the gas is given an additional velocity dispersion 
of Vvir/v^ in the spherical model, and V v i r /3 in the flattened, rotating model. This additional 
velocity dispersion is simply added quadratically to the thermal one, which would be valid if the 
gas were in optically thin clumps moving at this velocity dispersion. In practice, any relevant gas 
clumps in DLAs will be optically thick to Lya photons, but their inclusion would make our model 
much more complex. We will discuss the effect we would expect from optically thick clumps in §5. 

There are various sources to produce Lya photons: internal dissipation, star formation, and 
fluorescence caused by the intergalactic UV background. The external UV background will give rise 
to Lya fluorescence; at ~ 10 4 K, in an optically thick medium, recombination of the photoionized 
hydrogen has a 68% probability of producing a Lya photon (Spitzer 1978). Shock-heating of the 
gas will dissipate kinetic energy into heat, which will result in Lya emission after line excitation and 
collisional ionization followed by recombination. In addition, internal star formation can of course 
also produce Lya emission. We consider only fluorescent emission from the external background in 
this paper. In this case, the emissivity is proportional to the square of the ionized gas density. The 
other two sources of emission are more centrally concentrated, and therefore similar differences as 
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the ones between uniform emissivity and central source in the previous section can be expected. 

We take self-shielding into account in the spherical case, calculating the neutral hydrogen 
density of the cloud at every radius, where the total gas density profile is singular isothermal. 
We use an iterative algorithm similar to the one presented by Tajiri & Umemura (1998), as 
described in Zheng & Miralda-Escude (2002). Here, we assume a UV background intensity of 
10~ 22 ergss -1 cm -2 s" 1 Hz^ 1 sr _1 , constant at all frequencies between the H I Lyman limit, ul, and 
the He II Lyman limit, 4zv L . The intensity is set to zero at frequencies above \v^. This is a good ap- 
proximation since most higher frequency photons are likely to be absorbed by He II before reaching 
the hydrogen self-shielded zone. Even though in reality some very energetic photons exist, which 
can penetrate into the self-shielded region, their effect on the profile of neutral hydrogen is very 
small due to their low intensity. The profile of neutral hydrogen that is obtained is similar to that 
shown in Zheng & Miralda-Escude (2002). In the case of the flattened, rotating gas distribution, 
we simply take the same neutral hydrogen density profile as for the spherical case and flatten it 
with an axis ratio of 0.5 (the numerical code we have developed for computing the self-shielding 
correction applies only to spherically symmetric systems). This is of course not exact, although a 
much better approximation than neglecting self-shielding altogether. 



4.2. Results of Simulations 
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Fig. 5. — Surface brightness profiles of Lya emission for the spherically symmetric cloud (in mag- 
nitudes per unit solid angle, with an arbitrary normalization). Thick and thin lines are with and 
without velocity dispersion, respectively (see the text). Dashed line is for optically thin cloud (with 
an arbitrary vertical shift). 
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The surface brightness profile for the spherically symmetric, non-rotating cloud is plotted 
in Figure 5. The thick solid line shows the case of including the fluid velocity dispersion, a = 
Vvir/y^ = 60kms _1 , while the thin solid line is the result with only the thermal velocity dis- 
persion, (kT/mn) 1 / 2 = 12.8kms _1 . The dashed line is what is obtained when self-shielding is 
not included (it has been shifted vertically). The velocity dispersion has practically no effect on 
the surface brightness profile, because the change in spatial diffusion of the photons is very small. 
The emissivity due to the recombination has a "hole" in the inner part of the cloud because of 
self-shielding. Therefore, we see a "core" in the surface brightness profile (with a slight decline of 
surface brightness toward the center in the inner part due to geometric effects). The outer part of 
the profile follows the optically thin case. The core of the surface brightness profile is a signature 
that the emission is due to fluorescence of the external background, since other sources of emission 
are more centrally concentrated. 

Because the ionizing photons of the external background extend only up to Avl in frequency, 
essentially all the photoionizations and subsequent recombinations occur in the outer region of the 
gas cloud with a hydrogen column density of Nm < 10 19 cm" 2 , the inverse of the photoionization 
cross section at frequency 4^. The damped absorption wings are not yet important at this column 
density, therefore the photons escape the cloud after a single large jump in frequency when scattered 
by an atom in the high-velocity tail of the Maxwellian distribution, as discussed in §3 for the case 
of Lyman limit systems. 




Fig. 6. — Two-dimensional spectra of the spherically symmetric, non-rotating cloud. Left panel 
is for a thermal velocity dispersion of (kT/niH) 1 ^ 2 = 12.8 kms -1 only, right panel includes a fluid 
velocity dispersion of 60 kms -1 (see the text). The contours next to each other are separated by 1 
magnitude. 
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Figure 6 shows the 2-dimensional spectra of the Lya photons escaped from the non-rotating 
cloud. We assume that the slit is large enough to include the entire cloud. The double-peaked 
distribution of escaped Lya photons, with the intermediate frequency range of essentially zero flux, 
appears at every spatial position. 

Including an additional velocity dispersion has the effect of separating the two peaks. This is 
true in our model because we consider the case where there is no correlation between the velocity of 
two atoms that are spatially close, valid only when any gas clumps moving coherently are so small 
that they are optically thin to Lya photons. In other words, the additional velocity dispersion is 
effectively thermal. In the presence of coherent motions of optically thick regions, the intermediate 
frequency range with a highly reduced flux disappears, as discussed in §3 in the examples of an 
expanding or contracting cloud. 
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Fig. 7. — Two-dimensional spectra of the oblate rotating cloud viewed face-on. Left panel is for a 
thermal velocity dispersion of 12.8 km s _1 only, right panel includes an additional velocity dispersion 
of 35kms _1 (see the text). The contours next to each other are separated by 1 magnitude. 

Spectra of the flattened rotating cloud depend on the viewing angle and slit orientation. The 
face-on spectrum, plotted in Figure 7, is very similar to that of the spherical static case. Figure 8 
shows the edge-on spectra, with the slit perpendicular (top panels) and parallel (bottom panels) 
to the rotational axis. The spectrum is averaged along the spatial coordinate perpendicular to the 
slit; in other words, we are assuming that the entire cloud is included in the slit, but that the 
wavelength dispersion is wide enough that the size of the cloud does not introduce any smoothing 
of the spectrum. When the slit is perpendicular to the rotation axis, the rotation curve pattern 
is clearly seen. The velocity shift of the spectrum is about the same as the rotational velocity, 
V TO t = Kir a/2 /3 = 85kms _1 , at the edge of the system [note that the unit used in the horizontal 
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Fig. 8. — Two-dimensional spectra of the oblate rotating cloud viewed edge-on, with the slit 
perpendicular (top panels) and parallel (bottom panels) to the rotation axis. Left and right panels 
are for cases of different velocity dispersions as explained in Fig. 7. The contours next to each other 
are separated by 1 magnitude. 
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axis is Also = (2kT/mH) 1 ^ 2 ^o/c = 18.2kms _1 (fo/c)]. The variation of the shape of the spectrum 
along the slit is easily understood by thinking of an analogy to the expanding and contracting cloud 
in §3. The Lya photons are produced on the outer ionized layer of the system, near the radius r ss 
where the hydrogen becomes self-shielding. The line-of-sight velocity at this radius is proportional 
to V Iot r p /r ss , where r p is the projected radius on the slit position. Therefore, the central trough 
of the spectrum shifts linearly with the projected radius. Moreover, the variation of the velocity 
along the line of sight causes the peak of emission that is further from the mean velocity of the 
system to be enhanced, for similar reasons as in the case of the expanding and contracting clouds 
discussed in §3. 

For the case of the slit parallel to the rotational axis, the 2-dimensional spectrum displays a 
symmetric double-peaked pattern, which essentially results from averaging along the equator the 
spectrum with different projected velocity. This averaging results in a less sharp central trough and 
smoothing of the two peaks. In all cases, the velocity dispersion increases the distance between the 
two peaks in the photon distribution and broadens the width of each peak. 

5. Summary and Discussion 

We have developed a Monte Carlo code of Lya photon scattering, to construct simulated 
images and spectra. The code is fully three-dimensional and adaptable to a cloud with any given 
Lya emissivity, neutral hydrogen density distribution, and bulk velocity field. Several simple cases 
have been presented to show the images and spectra expected from a uniform gas distribution. We 
have applied the code to model DLA systems illuminated by the intergalactic UV background to 
study the spatial and frequency distribution of emitted Lya photons. Self-shielding produces a core 
in the Lya surface brightness profile. For a spherical, static cloud, Lya photons have a double- 
peaked distribution which is symmetric around the line center frequency. An oblate, rotating 
cloud shows the same double-peaked distribution with the rotational pattern of increasing mean 
wavelength with the position on the slit. Observations of the fluorescent emission from damped 
Lya systems with large telescopes (Hogan & Weymann 1987; Gould & Weinberg 1996) could reveal 
the presence of such systems and measure a velocity gradient, providing a direct measurement of 
the rotation rate. The Lya image would also reveal if the emission is due to fluorescence from 
the external radiation (in which case we expect a large core of the surface brightness at the radius 
where the hydrogen becomes self-shielded) or due to internal gas dissipation or star formation. 

The observations of metal lines in DLA systems show multiple absorption lines (e.g., Prochaska 
& Wolfe 1997, 1998), suggesting that the gas is clumpy. In this paper, we have assumed instead that 
the gas has a smooth distribution. The case of enhanced velocity dispersion we have considered is 
valid only if the gas is in very small clumps that are optically thin to Lya photons. Clumps that can 
give rise to the observed metal lines are highly optically thick. The effect of a clumpy gas distribution 
is not difficult to imagine. If there is only one clump along a line of sight, then this clump will 
essentially produce the same spectral shape in its emission line, shifted to its velocity. Thus, if the 
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emission line could be observed with sufficient angular resolution to resolve the individual clumps, 
the characteristic spectral shape of Figure 6 would appear at every position, although shifted to 
the clump velocity. If the clumps are not resolved, the emission line shape is obviously smoothed, 
just like in the case of an expanding or contracting cloud shown in Figure 2. Our numerical code 
can be applied in the future to the results of detailed hydrodynamic simulations of gaseous halos in 
the process of forming galaxies in cosmological realizations to predict other observable signatures 
that can test the mechanisms by which energy dissipation of gas leads to galaxy formation. 

A possible effect that we have not considered in this paper is the absorption of Lya photons 
by dust. Dust seems to be present in at least some DLAs as indicated by the reddening of the 
background quasars that have DLAs along the line of sight (Fall, Pei, & McMahon 1989; Pei, Fall, 
& Bechtold 1991) and the abundance ratios of iron-peak elements which follow a dust depletion 
pattern (Savage & Sembach 1996; Vladilo 1998). The resonant nature of the scattering of Lya 
photons by neutral hydrogen atoms increases the probability for them to be absorbed by dust 
grains. This is probably relatively more important when the Lya emission originates from a star- 
forming region (e.g., Meier & Terlevich 1981; Chariot & Fall 1993), where dust is likely to be 
abundant. If Lya photons of DLAs are caused by fluorescent emission as cases we consider in 
this paper, the absorption by dust would not be very severe due to the low column densities of 
the ionized outer layer (Chariot & Fall 1991). Since in the presence of dust the escape of Lya 
photons depends on many factors, such as the H I column density, the distribution of Lya source, 
the property of dust, and the topology of the medium (e.g., Ferland & Netzer 1979; Hummer & 
Kunasz 1980; Neufeld 1990, 1991; Chariot & Fall 1991), detailed modeling of DLAs is necessary for 
dust absorption to be taken into consideration. 

The authors thank the referee Joop Schaye for helpful comments. We also thank Xuelei Chen 
and David Weinberg for useful discussions. This work was supported in part by NSF grant NSF- 
0098515. 



APPENDIX 

To choose the velocity along the direction of the incident photon of the atom responsible for 
the scattering, we need to generate a random number u that has the following distribution (see 
eq.[4]): 

I{U) * {x-uf + a* (A1) 

where x and a are given. We give a brief description on how we generate random numbers with 
this distribution. Since the distribution has the symmetry under the transformation of (x, u) — > 
(— x, —u), we limit the following discussion to x > 0. 

We make use of the rejection method (Press et al. 1992) to generate u. The comparison 
distribution can be chosen to be g(u) oc [(x — u) 2 + a 2 ]" 1 , which can be integrated and inverted 
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analytically. A value of u with the distribution g(u) is first generated. We keep this value only if 

2 

a second random number uniformly distributed between and 1 is smaller than e~ u . In this way, 
we obtain values of n with the distribution f(u). In practice, when i > 1, the above method is 
inefficient because a very small fraction of values of u are not discarded. To increase the fraction 
of acceptance, we modify the comparison distribution to be 

/ \ I \{x — u) 2 + a?Y l , U < Itn , . „ N 

g(u) oc { ,, 2r / '„ , - (A2) 

w \ e~ u o[(x-u) 2 + U 2 ]' 1 , n>n v ' 

The value of no can be chosen to minimize the fraction of generated values that will be discarded. 

2 2 2 

The acceptance fractions are then required to be e~ u and e~ u /e~ u ° in the regions u < uq and 
u > no, respectively. 

A first random number R uniformly distributed between and 1 determines which region we 
use by comparing it with p, where 

P = = («. + \ )[(1 - + (1 + % = tan- 5!Z£. (A3) 

J_oo 9{u)du * 2 a 

We generate u through n = a tan + x, where 9 is a random number uniformly distributed in 
[— §,#o] and [#o, f ] for i? < p and R > p, respectively. Then another random number uniformly 
distributed between and 1 determines whether the generated value of n is accepted by comparing 
it with the corresponding fraction of acceptance. 
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